******************************************************************;
***figure8_table4.sas                                          ***;
******************************************************************;
***This program creates data for Figure 8 and Table 4 of the   ***;
***Davis, Grim, Haltiwanger, and Streitwieser 2013 RESTAT      ***;
***paper.                                                      ***;
******************************************************************;
***Output:                                                     ***;
*** ~/csv/figure8_pe.csv                                       ***;
*** ~/csv/figure8_tvs.csv                                      ***;
*** ~/csv/table4.csv                                           ***;
******************************************************************;

options ls=120 ps=max;

****MACROS****;

%macro v(var);

	s_&var = &var - mean_&var + gmean_&var;
	
%mend v;

%macro l(var);

data l&var;
	length vname $30.;
	set g&wt (keep = astd_&var gstd_&var gmean_&var);
	vname = "&var";
	rename astd_&var = astd_&wt;
	rename gstd_&var = gstd_&wt;
	rename gmean_&var = gmean_&wt;

%mend l;

%macro m(wt);

%*Across years;

proc summary data=r vardef=wdf;
	weight &wt;
	output out=g&wt
	mean(pf1) = gmean_pf1
	mean(dpf1) = gmean_dpf1
	mean(share_gen_hy) = gmean_share_gen_hy
	mean(share_gen_nu) = gmean_share_gen_nu
	mean(share_gen_pn) = gmean_share_gen_pn
	mean(share_gen_cl) = gmean_share_gen_cl
	mean(share_gen_ot) = gmean_share_gen_ot
	mean(puc_elect) = gmean_puc_elect
	mean(fac_i_adj) = gmean_fac_i_adj
	mean(int_pf1_fac) = gmean_int_pf1_fac
	mean(int_dpf1_fac) = gmean_int_dpf1_fac
	mean(rank_s_adj) = gmean_rank_s_adj
	mean(index1) = gmean_index1
	mean(index1_g) = gmean_index1_g
	mean(index1_yr) = gmean_index1_yr
	mean(index1_yr_pe) = gmean_index1_yr_pe
	mean(index1_yr_tvs) = gmean_index1_yr_tvs
	mean(mean_elas_nf_pe) = gmean_mean_elas_nf_pe
	mean(mean_elas_pe) = gmean_mean_elas_pe
	mean(mean_elas_nf_tvs) = gmean_mean_elas_nf_tvs
	mean(mean_elas_tvs) = gmean_mean_elas_tvs

	std(pf1) = gstd_pf1
	std(dpf1) = gstd_dpf1
	std(share_gen_hy) = gstd_share_gen_hy
	std(share_gen_nu) = gstd_share_gen_nu
	std(share_gen_pn) = gstd_share_gen_pn
	std(share_gen_cl) = gstd_share_gen_cl
	std(share_gen_ot) = gstd_share_gen_ot
	std(puc_elect) = gstd_puc_elect
	std(fac_i_adj) = gstd_fac_i_adj
	std(int_pf1_fac) = gstd_int_pf1_fac
	std(int_dpf1_fac) = gstd_int_dpf1_fac
	std(rank_s_adj) = gstd_rank_s_adj
	std(index1) = gstd_index1
	std(index1_g) = gstd_index1_g
	std(index1_yr) = gstd_index1_yr
	std(index1_yr_pe) = gstd_index1_yr_pe
	std(index1_yr_tvs) = gstd_index1_yr_tvs
	std(mean_elas_nf_pe) = gstd_mean_elas_nf_pe
	std(mean_elas_pe) = gstd_mean_elas_pe
	std(mean_elas_nf_tvs) = gstd_mean_elas_nf_tvs
	std(mean_elas_tvs) = gstd_mean_elas_tvs
	;

%*By year;

proc summary data=r vardef=wdf;
	by year;
	weight &wt;
	output out=m&wt
	mean(pf1) = mean_pf1
	mean(dpf1) = mean_dpf1
	mean(share_gen_hy) = mean_share_gen_hy
	mean(share_gen_nu) = mean_share_gen_nu
	mean(share_gen_pn) = mean_share_gen_pn
	mean(share_gen_cl) = mean_share_gen_cl
	mean(share_gen_ot) = mean_share_gen_ot
	mean(puc_elect) = mean_puc_elect
	mean(fac_i_adj) = mean_fac_i_adj
	mean(int_pf1_fac) = mean_int_pf1_fac
	mean(int_dpf1_fac) = mean_int_dpf1_fac
	mean(rank_s_adj) = mean_rank_s_adj
	mean(index1) = mean_index1
	mean(index1_g) = mean_index1_g
	mean(index1_yr) = mean_index1_yr
	mean(index1_yr_pe) = mean_index1_yr_pe
	mean(index1_yr_tvs) = mean_index1_yr_tvs
	mean(mean_elas_nf_pe) = mean_mean_elas_nf_pe
	mean(mean_elas_pe) = mean_mean_elas_pe
	mean(mean_elas_nf_tvs) = mean_mean_elas_nf_tvs
	mean(mean_elas_tvs) = mean_mean_elas_tvs

	std(pf1) = std_pf1
	std(dpf1) = std_dpf1
	std(share_gen_hy) = std_share_gen_hy
	std(share_gen_nu) = std_share_gen_nu
	std(share_gen_pn) = std_share_gen_pn
	std(share_gen_cl) = std_share_gen_cl
	std(share_gen_ot) = std_share_gen_ot
	std(puc_elect) = std_puc_elect
	std(fac_i_adj) = std_fac_i_adj
	std(int_pf1_fac) = std_int_pf1_fac
	std(int_dpf1_fac) = std_int_dpf1_fac
	std(rank_s_adj) = std_rank_s_adj
	std(index1) = std_index1
	std(index1_g) = std_index1_g
	std(index1_yr) = std_index1_yr
	std(index1_yr_pe) = std_index1_yr_pe
	std(index1_yr_tvs) = std_index1_yr_tvs
	std(mean_elas_nf_pe) = std_mean_elas_nf_pe
	std(mean_elas_pe) = std_mean_elas_pe
	std(mean_elas_nf_tvs) = std_mean_elas_nf_tvs
	std(mean_elas_tvs) = std_mean_elas_tvs

	median(pf1) = median_pf1
	median(dpf1) = median_dpf1
	median(share_gen_hy) = median_share_gen_hy
	median(share_gen_nu) = median_share_gen_nu
	median(share_gen_pn) = median_share_gen_pn
	median(share_gen_cl) = median_share_gen_cl
	median(share_gen_ot) = median_share_gen_ot
	median(puc_elect) = median_puc_elect
	median(fac_i_adj) = median_fac_i_adj
	median(int_pf1_fac) = median_int_pf1_fac
	median(int_dpf1_fac) = median_int_dpf1_fac
	median(rank_s_adj) = median_rank_s_adj
	median(index1) = median_index1
	median(index1_g) = median_index1_g
	median(index1_yr) = median_index1_yr
	median(index1_yr_pe) = median_index1_yr_pe
	median(index1_yr_tvs) = median_index1_yr_tvs
	median(mean_elas_nf_pe) = median_mean_elas_nf_pe
	median(mean_elas_pe) = median_mean_elas_pe
	median(mean_elas_nf_tvs) = median_mean_elas_nf_tvs
	median(mean_elas_tvs) = median_mean_elas_tvs

	p10(pf1) = p10_pf1
	p10(dpf1) = p10_dpf1
	p10(share_gen_hy) = p10_share_gen_hy
	p10(share_gen_nu) = p10_share_gen_nu
	p10(share_gen_pn) = p10_share_gen_pn
	p10(share_gen_cl) = p10_share_gen_cl
	p10(share_gen_ot) = p10_share_gen_ot
	p10(puc_elect) = p10_puc_elect
	p10(fac_i_adj) = p10_fac_i_adj
	p10(int_pf1_fac) = p10_int_pf1_fac
	p10(int_dpf1_fac) = p10_int_dpf1_fac
	p10(rank_s_adj) = p10_rank_s_adj
	p10(index1) = p10_index1
	p10(index1_g) = p10_index1_g
	p10(index1_yr) = p10_index1_yr
	p10(index1_yr_pe) = p10_index1_yr_pe
	p10(index1_yr_tvs) = p10_index1_yr_tvs
	p10(mean_elas_nf_pe) = p10_mean_elas_nf_pe
	p10(mean_elas_pe) = p10_mean_elas_pe
	p10(mean_elas_nf_tvs) = p10_mean_elas_nf_tvs
	p10(mean_elas_tvs) = p10_mean_elas_tvs

	p90(pf1) = p90_pf1
	p90(dpf1) = p90_dpf1
	p90(share_gen_hy) = p90_share_gen_hy
	p90(share_gen_nu) = p90_share_gen_nu
	p90(share_gen_pn) = p90_share_gen_pn
	p90(share_gen_cl) = p90_share_gen_cl
	p90(share_gen_ot) = p90_share_gen_ot
	p90(puc_elect) = p90_puc_elect
	p90(fac_i_adj) = p90_fac_i_adj
	p90(int_pf1_fac) = p90_int_pf1_fac
	p90(int_dpf1_fac) = p90_int_dpf1_fac
	p90(rank_s_adj) = p90_rank_s_adj
	p90(index1) = p90_index1
	p90(index1_g) = p90_index1_g
	p90(index1_yr) = p90_index1_yr
	p90(index1_yr_pe) = p90_index1_yr_pe
	p90(index1_yr_tvs) = p90_index1_yr_tvs
	p90(mean_elas_nf_pe) = p90_mean_elas_nf_pe
	p90(mean_elas_pe) = p90_mean_elas_pe
	p90(mean_elas_nf_tvs) = p90_mean_elas_nf_tvs
	p90(mean_elas_tvs) = p90_mean_elas_tvs
	;

%*Calculate standard deviation removing year effects and adding back in grand mean (year and grand mean unweighted);

data r&wt;
	merge r mnowt (keep = year mean_:);
	by year;

data r&wt;
	if _N_ = 1 then set gnowt;
	set r&wt;
	%*Subtract out annual mean and add in grand mean;
	%v(pf1)
	%v(dpf1)
	%v(share_gen_hy)
	%v(share_gen_nu)
	%v(share_gen_pn)
	%v(share_gen_cl)
	%v(share_gen_ot)
	%v(puc_elect)
	%v(fac_i_adj)
	%v(int_pf1_fac)
	%v(int_dpf1_fac)
	%v(rank_s_adj)
	%v(index1)
	%v(index1_g)
	%v(index1_yr)
	%v(index1_yr_pe)
	%v(index1_yr_tvs)
	%v(mean_elas_nf_pe)
	%v(mean_elas_pe)
	%v(mean_elas_nf_tvs)
	%v(mean_elas_tvs)

proc summary data=r&wt vardef=wdf;
	weight &wt;
	output out=astd&wt
	std(s_pf1) = astd_pf1
	std(s_dpf1) = astd_dpf1
	std(s_share_gen_hy) = astd_share_gen_hy
	std(s_share_gen_nu) = astd_share_gen_nu
	std(s_share_gen_pn) = astd_share_gen_pn
	std(s_share_gen_cl) = astd_share_gen_cl
	std(s_share_gen_ot) = astd_share_gen_ot
	std(s_puc_elect) = astd_puc_elect
	std(s_fac_i_adj) = astd_fac_i_adj
	std(s_int_pf1_fac) = astd_int_pf1_fac
	std(s_int_dpf1_fac) = astd_int_dpf1_fac
	std(s_rank_s_adj) = astd_rank_s_adj
	std(s_index1) = astd_index1
	std(s_index1_g) = astd_index1_g
	std(s_index1_yr) = astd_index1_yr
	std(s_index1_yr_pe) = astd_index1_yr_pe
	std(s_index1_yr_tvs) = astd_index1_yr_tvs
	std(mean_elas_nf_pe) = astd_mean_elas_nf_pe
	std(mean_elas_pe) = astd_mean_elas_pe
	std(mean_elas_nf_tvs) = astd_mean_elas_nf_tvs
	std(mean_elas_tvs) = astd_mean_elas_tvs
	;

data g&wt (keep = astd_: gmean_: gstd_:);
	if _N_ = 1 then set astd&wt;
	set g&wt;


%*Make long ds for output to csv;

%l(pf1)
%l(dpf1)
%l(share_gen_hy)
%l(share_gen_nu)
%l(share_gen_pn)
%l(share_gen_cl)
%l(share_gen_ot)
%l(puc_elect)
%l(fac_i_adj)
%l(int_pf1_fac)
%l(int_dpf1_fac)
%l(rank_s_adj)
%l(index1)
%l(index1_g)
%l(index1_yr)
%l(index1_yr_pe)
%l(index1_yr_tvs)
%l(mean_elas_nf_pe)
%l(mean_elas_pe)
%l(mean_elas_nf_tvs)
%l(mean_elas_tvs)

data l&wt;
	set lpf1 ldpf1 lshare_gen_hy lshare_gen_nu lshare_gen_pn
		lshare_gen_cl lshare_gen_ot lpuc_elect lfac_i_adj
		lint_pf1_fac lint_dpf1_fac lrank_s_adj lindex1
		lindex1_g lindex1_yr lindex1_yr_pe lindex1_yr_tvs
		lmean_elas_nf_pe lmean_elas_pe
		lmean_elas_nf_tvs lmean_elas_tvs
		;

proc sort data=l&wt;
	by vname;

proc print data=l&wt;
	title1 "Check l&wt";

%mend;

****START PROGRAM****;

data r;
	set elec.styr_elas; *Dataset created in figure7.sas;
        *Drop Nebraska - missing many regulatory variables;
        if postalst ne 'NE';

	int_pf1_fac = pf1 * fac_i_adj;
	int_dpf1_fac = dpf1 * fac_i_adj;

	nowt = 1;

proc sort data=r;
	by year;

*Create weighted summary stats;

%m(nowt)
%m(sum_pe)
%m(sum_tvs)

proc export data=msum_pe outfile="csv/figure8_pe.csv"
	dbms=csv replace;

proc export data=msum_tvs outfile="csv/figure8_tvs.csv"
	dbms=csv replace;


data l;
        merge lsum_pe lsum_tvs;
        by vname;

proc export data=l outfile="csv/table4.csv"
        dbms=csv replace;


run;
